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Abstract 


A new, definitive, reliable and fast iterative method is described for 
determining the geometrical properties of a shock (i. e. © Bn , n, V s and M^), 
the conservation constants and the self-consistent asymptotic magnetofluid 
variables, that uses the three dimensional magnetic field and plasma 
observations. The method is well conditioned and reliable at all 0g n angles 
regardless of the shock strength or geometry. Explicit proof of 
’uniqueness' of the shock geometry solution by either analytical or 
graphical methods is given. The method is applied to synthetic and real 
shocks, including a bow shock event and the results are then compared with 
those determined by preaveraging methods and other iterative schemes. A 
complete analysis of the confidence region and error bounds of the solution 


is also presented. 



1. Introduction 


The identification of an observed discontinuity as a shock rests on 
certifying a sequence of conditions (described below) which can only be 
rigorously expressed by specifying the geometrical orientation and speed of 
propagation of the discontinuity [Burlaga, 1971; Greenstadt et al . , 1984]. 
Within the various classes of shocks, there are diverse geometrical, 
theoretical and observational regimes which further differentiate shocks 
into quasi-perpendicular versus quasi-parallel, subcritical versus 
supercritical, laminar versus turbulent and resistive versus dispersive 
[Greenstadt et al., 1984; Edminston and Kennel, 1985; Kennel et al., 1985]. 
To specify which of these shock regimes a given set of observations 
illustrates is the initial task for the increasingly quantitative shock 
studies made possibly by the ISEE and Voyager instrumentation. Furthermore, 
the importance of the shock geometry in relation to the origin of particle 
reflection and acceleration of thermal particles near shocks has been 
emphasized by Sonnerup [1969], Gosling et al. [1982], Armstrong et al. 
[1985], Forman and Webb [1985], Wu [1984] and Goodrich and Scudder [1984] 
among others. Direct spacecraft measurements only determine these 
geometrical properties implicitly; they must be empirically inferred by 
solving what we shall describe below as the " Rankine-Hugoniot (RH) problem". 
This problem consists of taking the spacecraft observations (or time series) 
of density, velocity and magnetic field across a shock and finding a 
suitable Galilean frame where the discontinuity is time-stationary and where 
defensibly conserved quantities can be defined such as the normal mass flux, 
tangential stress, normal component of the magnetic field and tangential 
electric field together with the upstream and downstream asymptotic 
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magnetofluid states. The solution of the RH problem is a non-trivial 
problem involving specifying eleven, non-linearly intertwined, free 
variables. The angle of shock propagation Gg^ relative to the asymptotic 
upstream magnetic field § and the strength of the discontinuity 
characterized by the various Mach numbers can only be obtained after this 
frame shift, conservation constants and asymptotic states are determined. 
Once such a frame shift and states are found, a sequence of supplementary 
tests can be performed to determine if the discontinuity is a shock. In 
many respects it is easier to deny that a discontinuity is a shock than to 
guarantee that it is one. 

The single spacecraft determination of shock normals on planetary and 
interplanetary shocks have been previously studied by Colburn and Sonnett 
[1966], Chao [1970], Lepping and Argentiero [1971], Lepping [1972], 
Abraham-Shrauner [1972], Abraham-Shrauner and Yun [1976], Chao and Hsieh 
[1984] and Acuna and Lepping [1984] among others. Four basic methods of 
single spacecraft shock normal determination are widely used. These are 
magnetic coplanarity (MC), velocity coplanarity (VC), the least squares 
method of Lepping and Argentiero (LA) and the mixed data methods of 
Abraham-Shrauner (AS). All these methods use a subset of the RH 
conservation equations. The subset of the RH equations that restate 
conservation of mass flux, normal component of the magnetic field, 
tangential component of the momentum flux and tangential component of the 
electric field do not discriminate between four MHD discontinuous classes 
such as contact discontinuity, rotational discontinuity, tangential 
discontinuity and shock. Except for the Lepping and Argentiero method and 
the procedure described in this paper, the other approaches have used an 
even smaller subset of the above set to estimate the shock normal, speed and 
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geometry. These different methods have often revealed disparate results in 
the shock parameters estimated for the same set of observations. One of the 
difficulties on relying in some of these methods is that their use requires 
that one predefine the asymptotic magnetofluid variables by an "ad hoc" 
pre-averaging procedure. It is not clear in the presence of waves or random 
fluctuations that this kind of "ad hoc" procedure can describe the 
self-consistent asymptotic states of a shock. Alternatively, iterative 
schemes such as the LA method have tried to resolve this problem by solving 
directly for the asymptotic magnetofluid variables. These variables are 
subsequently used together with the magnetic coplanarity and mass flux 
conservation expressions to determine the shock normal and speed. Although, 
this approach is self-consistent, the LA method has the unfortunate 
difficulty that its 11-dimensional space of unknown magnetofluid variables ( 
p.|, pp, Vp - V 1 , and § 2 ), which span the parameter space, is large and 
irreducible. Besides, this 1 1 -dimensional space of variables is highly 
non-linear, giving rise for concern of the 'uniqueness' of the selected 
solution. Up to the present time, the problem of 'uniqueness' of the 
solution determined has remained completely unaddressed. Notice that 

methods that preaverage the data obviate questions of uniqueness by de facto 
algebraic computation. 

Another method used in the estimation of shock parameters has been the 
use of observations from two or more spacecraft [Chao, 1970; Ogilvie and 
Burlaga, 1969; Russell et al . , 1983a ,b]. Since situations where shock 

observations at more than one spacecraft are uncommon, we shall limit our 
discussion to comparisons with single spacecraft methods. In the situation 
where such observations are available we shall report the results of 
parameters determined from such multiple spacecraft methods. 
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This paper presents a new, fast iterative approach and solution to the 
RH problem to determine the shock parameters by means of a non-linear least 
squares method. An essential concept of this new method is that there exist 
a simple set of 'natural* variables that is separable. The new set of 
variables ( <j>, 6, V g , G n> B n , § t> , p 1 and p 2 ) together with two 

constraint equations form a vector basis that spans the 11 dimensional space 
of unknown parameters to be determined. By the term 'natural' we mean that 
choice of variables for which 'uniqueness' (or lack thereof) of the selected 
value is demonstrable either analytically (for linear variables) or 
graphically (for non-linear variables). Similarly, by separable we mean 
that the full set of 'natural' variables can be obtained through a 
self-consistent sequence of least squares problems each of which contains a 
small dimensional subspace (i. e. less than 2) of the complete set of 

unknown parameters. The existence of such an ordered sequence of smaller 
dimensional problems is a consequence of the fact that the RH equations 
which represent the model can be written in various forms permitting some of 
the unknown parameters to appear either explicitly or implicitly in the 

equations for the same set of observations. A further advantage of this 
approach over previous methods is that the nunber of linear variables of the 
unknown parameter space is large (i. e. seven) resulting in only four 
non-linear parameters of the full 1 1 dimensional space which require 

graphical 'uniqueness' investigation. By virtue of this separability we can 
explore the 'uniqueness' (or lack thereof) of all the possible minima that 
encompasses the optimal solutions of the RH equations. It is clear that the 
set of RH equations can support in addition to the shock solution other 
types of discontinuous solutions such as rotational, tangential and contact 
discontinuities which are inherent to the system of equations [Landau and 
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Lifshitz, I960; Burlaga, 1971; Akhiezer et al . , 1975]. After a thorough 
inspection of each of these minima and using a series of supplementary tests 
we can determine the most likely physical shock solution (if it exists) to 
the problem. Among the necessary conditions that an observed discontinuity 
must satisfy to be identified as a possible shock are: a)that in the 
selected Galilean frame, there should exist a defensibly non-vanishing mass 
flux, b)that there is a density and total electron plus ion temperature (if 
available) jump in the same sense across the discontinuity, c)that there 
should be a decrease of the normal component of the fluid velocity in the 
direction the density increases and c)that the predicted thermal normal 
pressure must increase with the density and should be comparable within the 
noise with the observed pressure (if available). 

In addition to providing explicit proofs of 'uniqueness' (or lack 
thereof) , the method converges equally fast for quasi-parallel or 
quasi-perpendicular shocks (for which extant methods converge extremely 
slowly requiring one/half day of VAX 11/780 computing time to determine a 
solution) . This new approach rarely takes more than a few seconds of 
computing time to correctly determine the shock parameters, the 
Rankine-Hugoniot conservation constants, as well as to graphically support 
the 'uniqueness' of the shock geometry selected. The method uses the 
observed plasma velocity and density as well as magnetic field measurements 
on both sides of the observed shock. The sequence of problems consists of 
initially determining the shock normal polar angles ( 4> , 9) and the shock 
speed V using the Rankine-Hugoniot equations and the plasma and magnetic 

3 

field data given by p, ^ and § on both sides of the shock by a non-linear 
least squares method. Once the optimal shock normal angles and speed have 
been determined , their values are used in conjunction with the data to 
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uniquely define the conserved constants. These constants are the mass flux 
G n , the normal component of the magnetic field B n , the tangential components 
of the momentum flux and the tangential components of the electric field 
in the frame of the observations. Finally, we use the determined shock 
normal, speed and conservation constants in conjuntion with the data back in 
the RH equations to predict the self-consistent asymptotic states of the 
magnetofluid in the upstream and downstream sides of the shock. We also 
estimate the error bounds and the region of confidence for the shock 
parameters. 

This paper is organized in the following manner. Section 2 presents a 
brief description of the RH conservation equations and their representation 
in an arbitrary reference system. The separable sequence of the least 
squares scheme for the solution of the shock geometry is presented in 
section 3- In section 4 we discuss the applications and results of this 
approach on simulated and real shocks. The results are then compared with 
those obtained by different techniques. Finally, a summary and conclusions 
of the results obtained is presented in section 5, with possible suggestions 
for future work. 

2. Rankine-Hugoniot Conservation Equations: The Model 

The determination scheme for the shock normal, shock speed, conservation 
constants and asymptotic states rests on a series of assumptions: 1)these 
parameters can be determined from the model equations of the 
Rankine-Hugoniot system; 2)there exist such a frame in which the shock is 
time-stationary; 3) the observations used in their determination constitute 
an ensemble of asymptotic states as predicted by the conservations 
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equations. This last assumption means that we are able to remove 
information associated with the shock (transition) layer. 

The conservations equations evaluated in the shock frame of reference 
(represented here by asterisk) for an isotropic plasma medium are [Boyd and 
Sanderson, 1969]: 


A [ pV * ] = 0 
n 

„ „ b i, 

A [ P V t - — — 1 = 0 

* 4 IT 

A[ B n (n x V t ) - V n (n x § t ) ] = 0 


(1) 

(2) 

(3) 


A[ B • n ] = 0 


(4) 


<B 2 - S>„ 2 > , 2 

A [ P + — + pV * ] = 0 

8ir 


(5) 


V * 2 

% 3£- 

a[ pV — + P v ( 


Y P B 


2 B ($* • 6) 


2 n (y-1 ) p 4irp 


) - 


] = 0 


(6) 


4 IT 


where p is the plasma mass density, is the plasma bulk velocity 

-> * 

component along the normal to the shock surface, is the flow velocity 

tangential to the shock surface, B and 6, are the associated normal and 

n l 

tangential components of the magnetic field, P is the total kinetic 

pressure, n is the normal unit vector and y is the adiabatic constant. The 

+ 

subscripts n and t imply projection operators defined for any vector A as A fi 
= I-n and = X- ( I-nn) where I is the unit tensor. The symbol A means that 

the quantity within brackets is to be evaluated after ('2') and before ( * 1 * ) 
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the shock transition layer as indicated by the time arrow and then 
substracted (i. e. AY = ^ ~ *■])• Equation (1) represents the mass flux 
conservation equation, (2) is the momentum flux conservation equation for 
the tangential components, (3) is the continuity equation for the tangential 
electric field, (4) is the continuity of the normal component of the 
magnetic field, (5) is the conservation of the normal momentum flux and 

•f -►# 

finally, (6) is the energy flux conservation equation. Note that B = B for 
a Galilean frame shift. If the plasma is anisotropic, equations (5) and (6) 
will change. In general the normal pressure term in (5) is represented by 
n*P*n where P is the full pressure tensor. However for an isotropic plasma, 
the tensor is diagonal and the expression n*P*n reduces to P. These system 
of equations can be simply expressed by means of a Galilean transformation 
into an arbitrary frame of reference (as for example the one where the 
observations are made) by the transformation 

t t (7) 

s 

where $ = V n represents the shock velocity and V is the plasma flow 

s s 

velocity in the frame of reference of the observations. 

It is clear from looking at these equations that they cannot be used 
without knowledge of the shock normal n, the shock speed V g as well as the 
quantities p, §, P and the constant y on each side of the shock. 

Equations (5) and (6) will not be used in our calculations. Although in 
some experiments the total kinetic pressure tensor (i. e. electron plus ion 
pressure) is known and in principle equation (5) could be used, this is not 
always so in all cases. Furthermore, in order for us to make a fair 
comparison of our method with others such as for example the LA method, we 
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shall restrict our system of conservation equations to (1) - (4) since they 

2 

are statements of proper conservation quantities within the approximation E 
2 

<< B . 

The Rankine-Hugoniot equations (1) - (4) writen in an arbitrary 
reference system using (7) are: 

A[G1 = A[p(tf - V, n) • n] = 0 (8) 

11 O 

A[B n ] = A[S • n] = 0 (9) 

(S*n) 

A[l] = A[p($ • n - V) (’&• ( I-nn) ) (§»(I-nn)) ] = 0 (10) 

S = 4ir 

A[l] = A[ n x ( V* ( I-nn) ) (S*n) - ($*n - V ) n x (S*(I-nn)) ] = 0 (11) 

b — 5 — 

where we have used V * = (^ - Vji)‘n and since V *(I-nn) vanishes 

in any arbitrary frame of reference. The variables G^, B n , and 
represent the conservation constants corresponding to mass flux, normal 
magnetic field, tangential momentum flux (stress) and tangential electric 
field, respectively. These equations represent a system of eight equations 
since (10) and (11) are vectorial expressions in an arbitrary system. In 
our notation the vector n = (n , n , n ) can be expressed in any orthonormal 

x y ^ 

system of coordinates where the observations are made, e. g. the 
heliocentric coordinate system (R, T, N). In addition to these equations we 
also have the normalization condition 

l"l = 1 

which acts as a constraint equation and allow us to reduce the space of 
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unknown parameters by one variable. This is accomplished by expressing the 
normal components in spherical coordinates as 

n^ = cose 

n^ = cos<|> sine where 0 < 4> < 2ir, and 0 < 6 < ir (12) 

n = sin* sine 
z 

Generally this selection of the sense of the shock normal direction is 
arbitrary. With these conditions the variables of the system are the two 

angles (<t>, 6), the shock speed V g and the magnetic field and plasma 
parameters. This final set of eight equations forms the basic system of 
equations that we use to determine the shock geometry by least squares. 

3. Application of the Sequence Method to the Rankine-Hugoniot Problem 

In this section we shall present the sequence of least squares problems 
that are used to determine the shock geometry using the model equations (8) 
- (11) described above. The basic RH problem can be stated as follows: 
given a typical ensemble of observations (i. e. a time series) of density p, 
velocity ^ and magnetic field 5 with random noise and/or waves superposed, 
characterizing disturbed states about a possible asymptotic (undisturbed) 
states and about a discontinuous change in fluid variables, estimate the 
optimal shock normal, shock speed, the conserved quantities across the shock 
and the appropriate compatible combination of magnetofluid variables that 
characterizes the self-consistent asymptotic states of the observed 
discontinuity. As initially posed, this problem requires the solution in an 
11 -dimensional space. In the subsequent sections it is shown that we can 
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reduce this multidimensional problem to a self-consistent sequence of least 
squares problems of smaller dimensions (i. e. less than or equal 2) each of 
which has a solution that can be demonstrated to be optimal, if not 
1 unique' . 

a. Shock normal and speed determination 

The first problem in the sequence is the calculation of the shock 
normal n and the shock speed V using equations (8) - (11) and the 
observations of density, velocity and magnetic field at both sides of the 
shock. We use a further simplification of the system of equations (8) - 
(11) by solving for the shock speed V g in equation (8) since it enters 
linearly in the equation. From equation (8) we get 

A[ptf] 

= • n (13) 

&[ p] 

which is the usual form for the shock speed. Substituting (13) into 
(9 )— (11) we get a system of seven equations in terms of density p, velocity 
magnetic field B and the shock normal polar angles (<t>, 0). If we now let 
the density, velocity and magnetic field to be given by the observations, 
then the final set of seven equations only contains two unknown parameters 
(<i>, 0). In our notation we define the model function ?'(x; p) = 0 as a 
vector of seven components formed by the final derived expressions. We also 
define the vectors x = (p^, p 2 , § 2 ) and p = (41, e) to represent 

the observations at both sides of the shock and the unknown parameters to be 
determined, respectively. Since P'(x; p) is a vector formed by the seven RH 
equations and we have a set of N pairs (i. e. both sides) of observations we 
now define a vector function ?(x; p) = { P\'(x..; p) } of size N' = 7 * N 
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which is the model function that represents all the observations since the 
index i varies fron one to N. The function F(x; p) can be expanded locally 
in Taylor series about some initial parameter set p^ = (<^°\ 6^) as 

? = m-, p (0) ) + — 
ap 

Equation (14) can be expressed in matrix notation as follows 

A Ap = AY (15) 

where a? = ? - ?(x; p^°^) is a vector of length N' where ? is the null 
vector (i. e. Y i = 0, i = 1, N') indicating that the conservation equations 
must be satisfied exactly. Equation (15) is called the normal equation of 
least squares. In this equation we also have defined A as a matrix of size 
N' x M (where M=2 the number of unknown parameters) formed by the partial 
derivatives of the seven model Rankine-Hugoniot equations with respect to 
the unknown parameters p = ($, 0) evaluated at the initial guess. These 
derivatives have been evaluated analytically and verified by numerical 
integration. 

For the sake of simplicity, we shall present the details of the least 
squares methodology of the solution of equation (15) in Appendix A. 
However, to summarize the results, the final solution of this equation is 
obtained by an iterative scheme that minimizes the norm of the residuals or 
the variance 0) = r^ r where r = (A Ap-AY). Once this minimum has 

been obtained the optimal solution p = p + Ap is recovered. From this 

+* » * + 

optimal solution p = ( <t> , 9 ) we can now recover the normal n using 

equation (12). We can easily demonstrate by graphical methods the 
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$ ft 

'uniqueness* of the solution p = (41 , 6 ) because there are only two free 

parameters. In this circumstance we can construct a plot of the contours of 
2 

the logx (e, 4>) function versus the angles e and <j> in the range 0 < e < ir 

and 0 < <f> < 2ir. Typical examples of the topology of these contour plots are 

displayed in Figures 2a-c for simulated shocks and Figures 5a-c for real 

shocks. These plots exhaustively illustrate the location in the polar shock 

2 

normal angles (4>, 6) where the logx (9. 4>) is a minimum. The gradient 
search selected location will correspond to a possible solution of the 
problem. However, if more than one minimum is present, each must be studied 
independently to determine the (<j>, 6) direction that is most consistent with 
the supplementary tests which characterizes a shock solution. If this 

situation occurs, additional information (as discussed in section 4) is 
required to ascertain the appr opiate physical shock solution. Once the 


normal n is obtained, equation (13) and the data are used to determine the 
optimal shock speed V . 

The solution for V s can be presented also as a one dimensional least 

squares problem. However, because V g enters linearly in equation (13) » it 

can be shown that the least squares problem has an analytic solution for the 

shock speed. To show this we write the least squares objective function 
2 

X (V ) that must be minimized 

s 


2 N A[p i V i ] - 2 2 

X CV ) = I ( - V ) 2 /o 2 

S i=1 AC P± ] 3 


where the first term in parentheses corresponds to the shock velocity as 
determined from the observations and a is the standard deviation of the 
shock speed obtained from the data. The value of n used in this calculation 
is the one obtained in the first part of this section. If we now take the 
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p 

first and second derivative of x (V ) with respect to V g it can be easily 
shown that the point V s where the minimum occurs is given by 

1 » *1.^1 „ 

V = Z • n 

N i = 1 A[p.] 

This is the only solution to the linear least squares problem in one 
dimension and its uniqueness can be analytically demonstrated by determining 
that the second derivative of x (V ) function is positive. The procedure in 

O 

Appendix A gives similar results since both approaches are mathematically 
equivalent . 


b. Determination of the Rankine-Hugoniot Constants 

We now proceed to determine the RH conservation constants. For this 
calculation it is convenient to rewrite the equations (8) - (11) in terms of 
the conservation constants as follows: 


= P( V • n - V s ) 

B = $ • n 
n 

(S*n) 

3. = p($ - V n) • n ($*(I-nn)) (S»(I-nn)) 

u S = i. = 


(16) 

(17) 

(18) 


c£. = n x (^* ( I-nn) ) (8 • n) - ($ • n - V ) n x (§*(I-nn)) (19) 

v — S — 


where the conservation constants G n , B n , and E t have been previously 
defined and c is the speed of light. An inspection of equations (16) - (19) 
indicates that these constants appear linearly and independently in the 
equations. This means that if we take the plasma and magnetic field in the 
right-hand side of equations (16) - (19) to be given by the measurements on 
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both sides of the shock, then the solution for the conservation constants 
reduces to a linear least squares problem whose solution can be obtained 
analytically* The method of obtaining these constants is similar to that 
previously used in the determination of the shock speed. Therefore the 
optimal conservation constants are given by 
1 2N 

G = I p (V * n - V ) (20) 

2N i=1 1 1 S 

1 2N 

B = — Z (S • n) (21) 

2N i = 1 

1 2N (S. -n) 

1 = I [ p. (1 • n - V ) (1 • (I-nn) ) - — (S .• (I-nn))] (22) 

t 2N i=1 1 1 3 i % 

1 2N 

ct = - — l [nx (^ . • (I-nn) ) (S. • n) - ($. • n - V ) nx(S • (I-nn))] (23) 
2N i=1 1 1 1 s i = 

These optimal solutions of the conservation constants are unique, 

2 

corresponding to an absolute minimum of the objective x function, since 
they result from a linear least squares problem whose second derivative can 
be shown to be positive. 

c. Determination of the self-consistent Rankine-Hugoniot asymptotic states 
In this section we proceed to determine the compatible RH asymptotic 
states. This is the final problem in the least squares sequence presented. 
Substitution of equations (16) and (17) into the vector equations (18) and 
(19) yield a set of six equations in terms of the conservation constants, 
the shock normal and the shock speed. After some algebraic manipulations, 


these six equations can be solved together to obtain the vector expressions 

15 



V(p) = 


(24) 


B 

G S + p — (n x ct ) 
C 4w L 


G ^ - p B *V(4ir ) 
n n 


+ n ( G /p + V ) 
n s 


p(B + G Cn x c&J) 

S< 0 ) 24 2-,= i * n B_ 


G n 2 - p B n 2 /(4.) 


(25) 


An inspection of these equations indicates that both the velocity and 
the magnetic field are functions only of the unknown density since the shock 
normal, speed and conservation constants have been previously determined. 
This implies that we only need to solve for the density at each side of the 
shock to predict the compatible Rankine-Hugoniot states. Two other 
important conditions that resulted naturally from (24) - (25) by taking the 
dot product of these equations with the shock normal n are 

n • = 0, n • = 0 (26) 


which means that in the frame of the observations, the product of the normal 
vector with the tangential momentum flux and electric field must vanish. 
This is not a surprising result since in the frame of the shock by 
definition these conditions must also be satisfied. These equations are 
singular for values of 


p = 0, 


4ttG_ 


The first condition (p = 0) represents an unphysical solution since for the 
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existanee of a shook the density at both sides must also exist. The second 
condition is more subtle and corresponds to solutions for which the 
asymptotic inflow speed is equal to the intermediate mode speed 

If* n-V 

v = 1-7 — f-l = 1 (2 ?> 

*A * n 

Here we have defined V A = B//(4irp ) as the Alfven velocity. This solution 

corresponds to a rotational discontinuity and not a shock. Notice that this 

is not inconsistent with the solution of the RH equations since a rotational 

discontinuity is also a solution to these equations. An inspection of the 

above equations shows that for any fast shock solution to exist the density 

2 2 

must lie in the range 0 < p < /B^ . In order for this regime to be 

physical requires that the mass flux G n should be experimentally non-zero. 

2 2 

For values of p > MirG^ VB^ the normal Alfven Mach number (M A ») is less than 
unity and this could indicate that either the disturbance is a slow shock or 
is not a shock at all. To assess whether the solution corresponds to a slow 
shock, additional information such as the temperature of both electrons and 
ions of the plasma is required. Another important consequence of this 
condition is that for perpendicular shocks (i. e. B n = 0) the singularity 
goes to infinity. This, of course, implies that only fast shock solution can 
exist in such conditions which is clearly compatible with MHD since slow 
shocks and rotational discontinuities becomes tangential discontinuities as 
0 Bn a PP roac h es 90° [Landau and Lifshiftz, I960]. 

To show the application of the least squares method to the system of 
equations (24) - (25) we again define a vector function P(x; p) representing 
the six equations given by (24)-(25) and one additional equation given by 
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the density observations as follows: 



( 28 ) 


The parameter = (p, 5)^ represents the plasma and magnetic field 

observations at either the upstream or downstream sides of the shock. We 
also define the index i which varies from 1 to the nunber of data points 
The variable p = (p) represents the unknown parameter to be obtained. As in 
the case of the shock normal angles, the function FC^; p) can be expanded 
in Taylor series (as in equation (14)) about an initial guess density value 
p^ = (p^ 0 ^) to give the expression 

A Ap = A? (29) 

which again represents the normal equation of least squares. In this case 
we define A as a matrix of size N' x 1 (where N' = 7 x N) formed by the 
partial derivatives of the seven model equations representing velocity V, 
magnetic field § and density p, evaluated at the initial guess. To avoid 
nunerical errors, these derivatives have been calculated analytically and 
verified by a nunerical quadrature. We also define AY as follows 
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-> 


AY 



- 5(p (0) )' 


S(p ( o) ) 


- P 


(o) 


representing the difference between the observations and the model 
parameters . 

The solution of the normal equation (29) is presented in Appendix A, 
however as in the shock normal angles situation, the final solution p = 

Hr 

(p ) of (29) is obtained by an iterative scheme that minimizes the variance 

2 -$-T •> 

X t p) = r r of the residuals. The determination of these asymptotic states 

can be divided in two parts. First, obtain by a least squares method the 

asymptotic state of the upstream side of the shock using the plasma and 

magnetic field observations, the previously determined shock normal, speed 

and the estimated conservation constants. By a similar procedure, the 

asymptotic states of the downstream side of the shock are determined. This 

approach is self-consistent since both sides of the shock yield the same 

conservation constants, shock normal and shock speed. To ascertain the 

'uniqueness 1 of the non-linear least squares iterative solution we can 

2 

simply graphically investigate the topology of the logx (p) function as a 

function of p at each side of the shock transition. In Figures 3a-c and 

Figures 6a-c we present examples of this function for simulated and real 

shocks, respectively. The solid and dashed lines represent the function 
2 # 

logx (p) versus density p (=p/p ) in the upstream and downstream sides, 
respectively. Because of the particular choice of the model equations (24) 
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2 2 2 

- (25) which are singular for values of p=0 and p = 4irG /I) , the v (p) 

n n * 

function has been pre-conditioned to discriminate against tangential, 

contact and rotational discontinuities. These RH solutions will correspond 

2 2 2 
to the maximim of the x (p) funtion for the singularity p = 4-irG^ /B^ . 
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4. Applications And Comparisons With Other Methods 


In this section we present the application and results of our method to 
both simulated and real shocks. We further compare the results obtained 
with those calculated by other techniques using the same data set. The 
simulated shocks were designed from the RH conservation equations [Tidman 
and Krall, 1971]. These shocks were constructed by prescribing the normal, 
the conservation constants, the shock speed and the angle between the 

shock normal and the upstream magnetic field. Once these parameters are 
specified the profiles of density, velocity and magnetic field were 
obtained. In an attempt to simulate the presence of waves or random noise 
in the data of an observed shock, the profiles of density, velocity and 
magnetic field were randomized independently. For simplicity, the random 
fluctuations superposed on these profiles were chosen to have a vanishing 
"time-average" with a relatively small amplitude (=10%). The final profiles 
were then used to evaluate and recover the shock parameters. We have 
selected a perpendicular (© Bn = 90°), parallel (© Bn = 0°) and oblique (© Bn = 
45°) synthetic shock as samples to test the method. We have also estimated 
the shock parameters for two real interplanetary shocks seen by the Voyager 
1 and 2 spacecrafts and a planetary bow shock crossing from the ISEE-1 
spacecraft. Comparison of our results with other methods including the two 
spacecraft method for the bowshock crossing are also presented. 

a. Synthetic Shocks 

Figures la-c show plots of the magnitude and components of the magnetic 
field and the plasma bulk velocity, together with the plasma density in an 
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arbitrary cartesian coordinate system of a perpendicular, parallel and 
oblique simulated shocks respectively. These shocks are designed to have a 

= 90°, 0° and 45° respectively with a shock speed of 500 km/sec. The 

4- 4>f 

perpendicular shock profile satisfies the condition B/p = constant and V = 
constant/p where the density profile is arbitrarily chosen to be 

p(t) = p + + p tanh[(t-t 0 )/x] 

where p + = (p 2 + p^/2 and p = (p 2 - p^/2, x controls the slope of the 
shock profile, t indicates the shock location and p^, p 2 are the asymptotic 
densities. The parallel shock velocity profile was chosen to be $ = 

constant/p and the magnetic field to be a constant across the transition 
zone. The density profile is chosen similar to the above expression for the 
perpendicular case. Although the synthetic shocks were constructed 
following Tidman and Krall [1971], we could have also designed them using 
equations (24) and (25) since they are equivalent. The oblique shock was 
designed following a new algorithm recently developed by Whang et al. [1985] 
which allow the construction of shocks for arbitrary 0g n angles (except 0° 
and 90°) given the plasma and magnetic field parameters in the upstream 
side . 

The vertical lines in Figures la-c indicates the data interval selected 
at both sides of the transition to evaluate the shock parameters. We now 
draw attention to the fact that there is no specific procedure on how and 
where the data should be selected. The only known requirement is that the 
data selected should not contain information about the transition layer, 
because in this region the RH conservation equations are not valid. 
However, there is no clear prescription on how far away from this layer or 
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how much data can be used to determine the shock geometry. Nonetheless, 
once the data interval , representing an ensemble of possible upstream and 
downstream states, has been decided; there is no restriction in either 
selecting equal or unequal number of data points at each side of the 
transition layer. Alternatively, since our method converges rapidly, we can 
select various data intervals with different number of data points to obtain 
an ensemble of solutions of the shock geometry. Then, we can investigate the 
intersection of all the solutions, within their error bounds, to 
statistically assess the shock geometry. 

The ’uniqueness® contour plot for the shock normal solution is presented 

for all the three cases in Figures 2a-c respectively. These figures show 

2 

the contour levels of the logarithm of the x objective function formed from 

all the data selected at both sides of the shock and the RH conservation 

equations, versus all the possible shock normal polar angles 0 and <|> as 

described in section 3a. Also indicated are shaded regions corresponding to 

2 

the lowest levels of the logx function, indicating the 95% confidence 
interval where the solution of the iterative scheme is located. Details 
about how to define such confidence intervals have been previously discussed 
by Scheffe [19593 and Bard [1974] and are presented in Appendix A. The 
topology of the 'uniqueness' surfaces of the perpendicular (Figure 2a) and 
the oblique shocks (Figure 2c) seems to be similar. Both surfaces show the 
solution to lie inside a "ridge" where the value of the contour levels are 
the smallest. However the topology of the parallel shock (Figure 2b) not 
only indicates the presence of a "ridge" but it shows a pair of "holes" at 
conjugate (i. e. anti-collinear) angles. It is important to note that the 
"holes" and "hills" shown in these topologies always appear in conjugate 
pairs due to the sign ambiguity in the shock normal solution. These 
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topologies seems to be typical of the type of shock in study, however this 
should be substantiated by a statistical study. 

A search for a solution through all the holes shown in these figures 
indicates that not all of them correspond to the proper shock solution of 
the problem. To assess the proper shock-like solution, four conditions must 
be considered. First, we should certify a defensibly non-vanishing mass 
flux G n (i. e. | 6 G n /G n | < 1). Secondly, we must compare the quality of the 
asymptotic magnetofluid states predicted with the corresponding observed 
variables and determine whether such predictions are within the standard 
deviations of the measurements. Thirdly, using the asymptotic magnetofluid 
variables we determine the Alfven Mach nunber (M^ = M^' cose^). If the 
quality of the asymptotic states is acceptable and the diagnosis of the 
problem indicates a fast shock solution, then the normal Alfven Mach nunber 
must be theoretically greater than unity. However, if the normal Alfven 
Mach number is computationally smaller than unity, suggesting the 
possibility of a slow shock, we must consider the relative mass flux error 
| 6 G n /G n | and the additional temperature information to correctly assess the 
final solution. Finally, using the RH equation (5) we may predict the 
thermal pressure junp across the shock given Dy 


A(n* P»n)= aP = — [ (— + p(tf«n-V ) 2 ) . - (— + p (tf.n-V ) 2 ) n ] 

8 * 8 * 

where the subscripts "d" and "u" represent the downstream and upstream sides 

respectively. Note that the prediction of the scalar normal pressure junp 

across the shock is independent of an assunption of an equation of state. 

To evaluate AP we use the predicted asymptotic magnetofluid variables. If 
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AP yield a negative pressure value, then the "hole" selected cannot 


correspond to a shock. 

The shock normal solution obtained by the pre-averaged and iterative 
schemes are also presented in Figures 2a-c. These solutions are indicated 
by various symbols corresponding to the method indicated in Tables la-c. In 
situations where different methods yield the same solution or very near each 
other, the symbol indicator corresponding to each technique will point with 
an arrow to the proper location in the contour plots to avoid overcrowding 
the solutions. 

To determine either the Bg^ angle defined by B^ = cos - ^ (lS Ue n/|^ u | ) , the 

Alfven Mach number (= M^' cosGg n ) , or the pressure jump condition, it is 

necessary to evaluate the asymptotic states. By evaluating the optimal 

density state at each side of the shock, the self consistent asymptotic 

velocity and magnetic field are 'uniquely' determined. Figures 3a-c show 

the 'uniqueness' of the solution for the evaluation of the asymptotic states 

2 

in all three cases. These figures indicate the levels of the logx (p) 

objective function formed from the data selected at both sides of the shock 

and the model equation (28) versus the normalized density p = p/p as 

described in section 3c. In Figures 3a-c the solid and dashed curves 

2 

represents the logx (p) f° r the upstream and downstream side of the 


transition layer, respectively. The normalization density p corresponds to 

the final predicted value determined by the iteration scheme at each side of 

the shock. For the three types of simulated shocks presented, these curves 

only contain a single minimum corresponding to the value correctly 

determined by the iteration scheme. The fact that only one minimum exists 

* 2 2 

at a density value p in the range 0 < p < %G n /B indicate, not only the 
'uniqueness' of the shock-like solution for the density, but also for the 
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asymptotic velocity $(p ) and magnetic field fi(p ) as described in equations 

(24) and (25). Note also that Figures 3a-c show the singular behavior of 
2 

the x (p) function when p = 0. We previously have established the existance 

2 2 

of another singularity p = 4irG n /B^ in section 3c which indicates the 

transition from fast to slow shock. This singularity also exists in these 

2 * 

cases, but they are located far away from the x (p ) minimum and off the 
figures. 

The general results for the three synthetic shocks are summarized in 
Tables la-c. These tables contain the results obtained from the 
pre-averaged and iterative methods for the geometrical characteristics of 
the shocks. For comparison, the first column contains the exact solution of 
the shock geometry of the synthetic shocks and the last column indicated by 
VS shows the solution obtained by our method. The first nine rows show the 
geometrical parameters that describe the shock geometry. The next fourteen 
rows show the asymptotic magnetofluid variables used by the pre-average 
methods and those determined from the iterative schemes. Finally, the last 
two rows give a measure of the efficiency of the iterative schemes in 
obtaining a solution. 

For the perpendicular shock (Figure 2a) the solution of our iterative 

scheme gives a 0g n = 90.0° ± 2° and a shock speed of about 503 ± 17 km/ sec. 

The final solution is located at 0 = 20° ± 0.1° and <t> = 160° ± 14.7° as 

indicated by the dark circle inside the shaded region (95% confidence 

interval) along a ridge in Figure 2a. At this location the final value of 
2 

the logx is -0.33. The path followed by the descending iterative gradient 
scheme has been indicated by the connected circles. Because of the sign 
ambiguity in the shock normal, a second solution exists at conjugate angles 
0 = 160° ± 0.1° and <(> = 340° ± 14.7°. This second solution represents the 
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normal vector opposite ( anti-collinear) to the one indicated in Table la and 
has a V g of opposite sign to its mirror image. Both solutions are perfectly 
valid, however in general the proper sign of the solution is decided by 
compensating the sign of the scalar shock speed V s with the obtained normal 
to form the vector shock velocity ^ = V n. Besides our solution we also 

show in Figure 2a the solution obtained by other methods. The general 
results of the analysis for this perpendicular shock are presented in Table 
la. A comparison with other methods of the results suggests that except for 
the LA method whose convergence to a solution was extremely slow and for the 
MC method, which did not reproduce the known solution accurately enough, all 
other procedures yield reasonable results relative to the exact solution. 
The reason the MC method gave a poor solution seems to be related to the 
fact that the MC equation 

+ (§ 1 xB 2 )x(B 1 -B 2 ) 

n = ± 

|(b 1 xb 2 )x(b 1 -b 2 )| 

becomes singular for perpendicular shocks. The convergence in the LA method 
was too slow because at each step in the iteration process it depends on the 
same expression of magnetic coplanarity (MC) [Lepping and Argentiero, 1971]- 
Besides, we noticed that the LA method is allowed to search for solutions in 
unphysical regions where, for example, the density is predicted negative. 
In consequence, this kind of unconstrained scheme slows down the iteration 
process and permits the gradient search to take large steps that may well 
violate the initial local Taylor series expansion. Recently Acuna and 
Lepping [1984] attempted to speed up the convergence rate of the LA method. 
Although sane increase in the rate of convergence was obtained, this has not 

controlled and constrained the search in un physical regimes. One important 
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aspect which arises from the results in Table la is that both the AS method 


given by 


_ [(S 2 -6 1 )x(^ 2 4 1 )]x(fe 1 -6 2 ) 

n = ± — 

|[(b 2 -b 1 )x(v 2 -v 1 )]x(b 1 -b 2 )| 

and the VC method given by 


kVV 1 

yield accurate solutions for the perpendicular shock geometry. Although the 
VC method is an approximate technique, in the case of high Mach number 
perpendicular shocks, it is theoretically expected to produce the proper 
solution as argued by Abraham-Shrauner and Yun [1976]. 

Another aspect which resulted from the analysis of the contours in 
Figure 2a is the existence of two unphysical minima located at conjugate 
pair of angles e = 90° ± 10° and $ = 70° and 250° ± 10°. This class of 
minima are almost always present in the contours. Their location are nearly 
orthogonal to the proper optimal solution. These solutions yields M^' < 1 
and the magnetofluid variables determined from them are in very poor 
agreement with the plasma and magnetic field observations. Moreover, these 
solutions violate the pressure junp condition across the shock layer. An 
inspection of these unphysical shock solutions suggests that they seem to 
belong to either the family of the tangential, contact or rotational 
discontinuities since the mass flux is numerically extremely small. 

A similar analysis has been performed for a synthetic parallel shock as 
shown in Figure 1b. The plots in Figures 2b and 3b show the 'uniqueness' of 
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the shock normal solution and the magnetofluid variables for this shock. 

Superposed on Figure 2b we show the locus of the descending path of the 

iterative scheme and the solutions by other methods. The general results of 

the analysis of these shocks are tabulated in Table 1b. The results of our 

iterative scheme yield 0^ = 0.03° ± 5° with a shock speed of 500 ± 12 

km/sec. The final solution of the parallel shock is indicated by the dark 

circle in Figure 2b at the polar angles e = 20° ± 11.8° and <j> = 160° ± 3k. 6° 

located in one of the isolated shaded "holes". At this location the final 
2 

value of logy is 1.28. Similarly there is a conjugate solution at e = 

160.1° ± 11.8° and $ = 3k0° ± 3k. 6° corresponding to the opposite sign of 

the shock normal. Both solutions are physical since they predict a positive 

pressure jump across the layer. For this case the LA method did not 

converge within a reasonable time. Furthermore, neither the MC method nor 

the AS method predicted the correct solution for this shock because of the 

singular behavior that both methods have as the Gg n approaches 0°. The 

location of the MC and AS solutions, shown in Figure 2b, indicates that they 

2 

reside in a deep shaded 'ridge', where the x function is small and where 
the unphysical solution exists. On the other hand the VC method gave good 
results that lies within the 95% confidence level of the minimum. However 
this agreement may be fortuitous because in the design of the parallel 
shock, the normal was chosen to be along the direction of the flow velocity 
which is a basic assumption of the VC method. 

The analysis of the oblique shock in Figure 1c yield a solution of e Bn = 
k5.2° ± 8° and a shock speed of k99.8 ± 19 km/sec. This final solution is 
indicated in Figure 2c by the dark circle along the shaded "ridge" at the 
polar angles 0 = 20.8° ± 1.0° and 4> = 159.3° * k6.0°. A conjugate solution 
is also present at 6 = 159.2° ± 1° and <j> = 339.3° ± k6.0°. The results of 
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PLASMA AND MAGNETIC FIELD FOR SYNTHETIC PERPENDICULAR SHOCK 
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Figure la. Magnetic field and plasma data plots of a synthetic perpendicular shock. Vertical 
lines indicates the data interval selected for the shock geometry analysis at 
both sides of the layer. The horizontal axis represents time in arbitrary units 
and the shock time is 0.5. 
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PLASMA AND MAGNETIC FIELD FOR SYNTHETIC PARALLEL SHOCK 
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Figure lb. Magnetic field and plasma data plots of a synthetic parallel shock. Vertical lines 
indicates the data interval selected for the shock geometry analysis at both 
sides of the layer. The horizontal axis represents time in arbitrary units and 
the shock time is 0.5. 
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Figure lc. Magnetic field and plasma data plots of a synthetic oblique shock. Vertical 

lines indicates the data interval selected for the shock geometry analysis at both 
sides of the layer. The horizontal axis represents time in arbitrary units and 
the shock time is 0.5. 
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Figure 2a. Contour plots of the log x 2 (0, $) function versus the shock normal polar angles 
( 9 , 0) indicating the ‘uniqueness’ of the shock geometry for a synthetic perpen- 
dicular shock. Superimposed on these figures we indicate the location of the so- 
lution by magnetic coplanarity (MC) ‘*\ velocity coplanarity (VC) *+’, Abraham- 
Shrauner (AS) ‘A’, Lepping-Argentiero (LA) ‘CT and our solution (VS) ‘©’. 
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Q (deg . ) from X 

Figure 2b. Contour plots of the log x 2 (6, <t>) function versus the shock normal polar angles 
(6, 0) indicating the ‘uniqueness’ of the shock geometry for a synthetic parallel 
shock. Superimposed on these figures we indicate the location of the solution by 
magnetic coplanarity (MC) ‘*\ velocity coplanarity (VC) ‘+\ Abraham-Shrauner 
(AS) ‘A’, Lepping-Argentiero (LA) and our solution (VS) *•’. 
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0(deg.) from X 


Figure 2c. Contour plots of the log x 2 (0, 0) function versus the shock normal polar angles 
(d, 0) indicating the ‘uniqueness’ of the shock geometry for a synthetic oblique 
shock. Superimposed on these figures we indicate the location of the solution 
by a magnetic coplanarity (MC) **’, velocity coplanarity (VC) *+’, Abraham- 
Shrauner (AS) ‘A’, Lepping-Argentiero (LA) ‘d’ and our solution (VS) ‘®\ 
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Figure 3 
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Plots of the log x 2 (P) function versus the normalized density p = p/p* indica- 
ting the ‘uniqueness’ of the asymptotic magnetofluid variables in the upstream 
and downstream sides of a synthetic (a) perpendicular, (b) parallel and (c) ob- 
lique shocks. The normalization constant p* is the value obtained by the 
iteration scheme at the minimum of the X 2 function for each side of the shock. 
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Table la. Results of the analysis of the synthetic perpendicular shock. 
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Table lb. Results of the analysis of the synthetic parallel shock. 
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Table lc. Results of the analysis of the synthetic oblique shock. 
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the analysis for this shock are summarized in Table 1c. From the results in 
Table 1c we note that the MC and AS methods yield reasonable solutions, 
however the VC solution resulted in poor agreement with the exact solution. 
This is due to the mis-alignment of the bulk velocity to the shock normal 
and also probably to the small value of the Alfven Mach number. The angle 
0g n and the shock speed V s given by the VC method are well outside the 
confidence bounds of the proper minimum solution. On the other hand, the LA 
method yielded good 6g n angle within the 95% confidence region; however, 
both the LA shock speed and pressure jump across the shock depart 
considerably from the exact solution. Note that the LA method cannot 
predict the asymptotic plasma bulk velocity at each side of the shock, but 
it can only resolve the velocity jump A^ = ^ - V-| across the layer. The 
predicted velocity jump across the layer obtained by the LA method for this 
case yielded AV = (115.8, -68.9, 21.6) km/sec which compares relatively well 
with the exact velocity jump AV = (111.4, -20.5, 7.4) km/sec obtained from 
Table 1c. Moreover, we find a conjugate pair of unphysical solutions at e = 
90° ± 5° and 4 > = 70° and 250° ± 5° that violate the pressure jump condition 
across the shock and are located almost orthogonal to the proper solution. 
This unphysical "holes" also give a very small mass flux suggesting that 
this candidate solution is either a tangential or a contact discontinuity. 

b. Real Interplanetary and Planetary Shocks 

Complete plots of the magnitude and components of the magnetic field and 
plasma bulk velocity, with the plasma density in a heliocentric (R, T, N) 
coordinate system for a quasi-perpendicular and a quasi-parallel 
interplanetary shock are presented in Figures 4a-b, respectively. A similar 
plot for a planetary quasi-perpendicular bow shock in a GSE (geocentric 
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solar ecliptic) coordinate system is also shown in Figure 4c. The intervals 
selected for the analysis of the shock geometry are indicated by the 
vertical lines. 

The event on November 27, 1977 corresponds to a quasi-perpendicular 
forward shock at 2225:57 UT as seen by the Voyager 1 at about 1.6 AU. 
Figures 5a and 6a represents the ’uniqueness' plots for the shock normal and 
asymptotic magnetofluid variables solutions respectively. Superposed on 
Figure 5a are the solutions obtained by other techniques. The path followed 
by the iterative gradient scheme to get to our solution is indicated by the 

connected circles. The results of this event are tabulated in Table 2a. 

First, note the similarity of the topology of the shock normal conjugate 
pair of solutions to that of the simulated perpendicular shock. As before, 
the solutions are along a ’’ridge" path and are located at two "thin" shaded 
contours centered about 0 = 37.5° ± 1.8° and $ = 262.5° ± 22° and 0 = 142.5° 
± 1.8° and <j> = 82.5° ± 22° where the value of logy 2 is -0.014. Our 

estimates confirm that this event is a quasi-perpendicular shock with a 0g n 
= 84.2° ± 9° and a shock speed of 305.5 ± 19 km/ sec. Comparison of our 

solution with those obtained by other methods is shown in Table 2a. An 

inspection of the asymptotic magnetofluid variables predicted by our method, 
compared to the average values used by the pre-averaged techniques, and 
their standard deviations as shown in the first column, indicates the good 
agreement of our predictions within the error bounds of the data. For this 
event the LA solution is unknown beacuse the method did not converge within 
a reasonable time. Nonetheless, both the AS and the VC methods yielded good 
solutions because the shock meets the preconditions of these methods. Both 
solutions lie within the 95$ confidence region about the minimum and they 
are within the error bounds supported by the data and the calculations. 
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However, although the MC solution is not extremely different from those 
obtained by other methods, it is nevertheless, outside the accepted 
confidence level. We have also estimated for this event using the electron 
and proton data, the observed thermal scalar pressure jump across the shock 
layer. The average electron temperature in the upstream and downstream 
sides of the shock are 6.0 ev and 11.0 ev respectively. Similarly, the 
proton temperatures in the upstream and downstream sides are 0.8 ev and 3-5 
ev, respectively. Assuning charge neutrality we find that the thermal 
pressure junp is about 224 ev/cc. The value predicted by our method (see 
Table 2a) gives 296 ev/cc. This discrepancy of about 30% in the prediction 
of A(n*P*n) can be explained by taking into consideration the geometry and 
orientation of the electron detector in the Voyager 1 spacecraft. The fact 
that there is only one electron detector which points always perpendicular 
to the radial direction almost in the equatorial plane (i. e. T-N plane) 
certainly indicates that the temperature reported are underestimated since 
there is not enough directional coverage of the electron distribution 
function to determine the proper pressure tensor. Besides, the important 
temperature component required for the pressure junp calculation should be 
that along the normal. But since this event is a quasi-perpendicular shock, 
this indicates that we must evaluate T A with certainty. An inspection of 
the electron detector orientation seems to indicate that the temperature 
obtained from it is the parallel component because of the field geometry 
relative to the detector during this period. 

The case on January 29, 1978 is a quasi-parallel reverse shock at 

0918:39 UT seen by Voyager 2 at about 2 AU. This shock has been previously 
studied in association with its structure by Scudder et al. [1984] and in 
the context of upstream waves by Vinas et al . [1984], The 'uniqueness" 
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plots for the shock normal and magnetofluid variables solutions resulting 

from our calculations are shown in Figures 5b and 6b respectively. The 

results by our technique and by the other methods are shown in Table 2b. 

For reference, the shock normal results of other methods are superposed in 

Figure 5b. We find and confirm that this event is a quasi-parallel shock 

with a 0g n = 29° ± 18.0° and a speed of 261 ± 39 km/sec. The shock normal 

corresponding to this event is located at 0 = 157.7° ± 14.1° and <j> = 125.9° 

± 31.1° with its conjugate normal at 0 = 22.3° ± 14.1° and <j> = 305.9° ± 

2 

31.1° where the value of the minimum logx is -0.9. Besides our solution, 
the AS method gives the only other result which lies within the error bounds 
of the accepted solution. All the other methods lie outside the 95% 
confidence interval. Note that the MC and LA methods are well outside the 
region where the minimum is located indicating that their solution are 
poorly resolved. Another important aspect of our calculations is the good 
agreement of the predicted thermal scalar pressure jump across the layer 
with the observed thermal pressure jump as obtained from the electron and 
proton data. The average electron temperature in the upstream and 
downstream sides of the shock are 6.1 ev and 6.2 ev respectively. The mean 
proton temperature in the upstream and downstream sides are 1.92 ev and 5.0 
ev respectively. Assuming equal density for electrons and ions 
(quasineutrality) we find that the thermal pressure jump is about 6.M ev/cc. 
Comparing this value with our prediction in Table 2b we find agreement well 
within the 10% error of the observed jump while that obtained by other 
methods are larger . 

The final event we investigated is a planetary bow shock crossing from 
ISEE-1 spacecraft on November 7, 1977. The shock crossing time is at 

2251:19 UT and the data intervals selected at each side of the layer for the 
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analysis are indicated by the vertical lines in Figure 4c. So far, we have 
used only proton plasma data to analyze the shock geometry. However, for 
this event we shall use the electron plasma data obtained from the Goddard 
three dimensional electron spectrometer . The 'uniqueness' plots for the 
shock normal and the asymptotic magnetofluid variables are shown in Figures 
5c and 6c respectively. As usual, the locus of the iterative scheme and the 
results from all the methods are indicated in Figure 5c. The overall 
results of the analysis of this event are presented in Table 2c. Our 
analysis indicates that this event is a quasi-perpendicular shock with 8^ = 
74.4° ± 20° with a shock speed of -8.4 ± 31 km/s. The solution is located 
inside one of the shaded "holes" representing the 95% confidence region at 
the polar angles 8 = 164.9° ± 7.5° and $ = 332.9° ± 28° where the minimum of 
the logx is 1.7. Another conjugate solution is also found at 0 = 15.1° ± 
7.5° and <f> = 152.9° ± 28° corresponding to the opposite normal sign selected 
in Table 2c. The solutions obtained by the AS and the VC methods are also 
very near the optimal minimum solution. Despite the fact that the AS and VC 
solutions are within the 95% confidence region, their relative shock speed 
error is greater than 10% compared to the shock speed determined from the 
two spacecraft method. However, the MC and LA methods yield very poor 
solutions, which are well outside the acceptable confidence interval. 
Indeed, the MC solution is quite close to one of the unphysical solutions of 
the problem. For comparison, the velocity jump across the shock determined 
by the LA method gives A V = (15.4, -56.3. 56.2) km/sec while our solution 
(VS) gives AV = (197.2, -46.0, 26.5) km/sec. Figure 5c also show the 

presence of a conjugate pair of unphysical shock solutions that yield 
negative pressure jump across the shock layer. These unphysical shock 
'solutions' are located at 8 = 76 ° ± 8° and <t> = 353° ± 9.5° and also at its 
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conjugate position e = 104° ± 8° and 4 > = 173° ± 9.5°. This bow shock has 

been exhaustively investigated by Scudder et al. [1985]. They have reported 

two spacecraft calculations of the shock speed using the ISEE-1 and -2 

observations of the same shock crossing. We have compared our calculations 

of the shock speed with that determined by the two spacecraft time delay 

method and the result is in excellent agreement with it. From the 

*► 

separation distance between the spacecrafts AS = (115.2, -193.0, 111.4) km, 
the time delay of the bow shock crossing At = 26 sec and assuming the shock 
normal determined by our method we can find the shock speed as seen by an 
observer in the spacecraft frame 


( spc) 


Therefore the shock speed y gives -8.8 km/sec. A comparison with our 

results indicates an excellent agreement within the error bounds of the 
calculations. Scudder et al . [1985] have also reported the velocity using a 
somewhat larger data interval in the downstream side of the shock. Their 
solution is also consistent within their error bounds with that determined 
in this paper . 

We have also evaluated and compared the thermal pressure jump A(n*£*n) 
across the shock with that calculated from the electron and proton 
temperature data for the data interval indicated in Figure 4c. The average 
electron temperature in the upstream and downstream sides of the shock are 
1.39 ev and 4.0 ev , respectively. Similarly, the proton temperatures in the 
upstream and downstream sides are 6.0 ev and 148.2 ev respectively. 
Assuming again charge neutrality we find that the observed thermal pressure 
jimp across the shock is 4736,1 ev/cc. The predicted pressure jump (see 
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VI -1977 MAGNETIC FIELD AND PLASMA DATA FOR SHOCK, NOV. 27, 1977 
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Figure 4a. Plasma and magnetic field data time plots for a real quasi-perpendicular inter- 
planetary shock seen by the Voyager 1 spacecraft. The vertical lines represent 
the data interval selected for the shock geometry analysis. 
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V2-1978 MAGNETIC FIELD AND PLASMA DATA FOR SHOCK, JAN. 29,1978 
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Figure 4b. Plasma and magnetic field data time plots for a real quasi-parallel interplane- 
tary shock seen by the Voyager 2 spacecraft. The vertical lines represent the 
data interval selected for the shock geometry analysis. 
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ISEE-1 PLASMA and MAGNETIC FIELD for BOWSHOCK NOV 7, 1977 



Figure 4c. Plasma and magnetic field data time plots for a real planetary bow shock seen 
by the ISEE-1 spacecraft. The vertical lines represent the data interval selected 
for the shock geometry analysis. 
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Figure 5a. ‘Uniqueness’ contour plots of the log x 2 ( 0 , 0) function versus the shock normal 
polar angles (6, 0) of a real quasi-perpendicular interplanetary shock. The loca- 
tion of the solution of magnetic coplanarity (MC) velocity coplanarity (VC) 
‘+\ Abraham-Shrauner (AS) ‘A’, Lepping-Argentiero (LA) *□’ and our solution 
(VS) ‘®’ are indicated. 
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UNIQUENESS OF NORMAL 



Figure 5b. ‘Uniqueness’ contour plots of the log x 2 (0, <P ) function versus the shock normal 
polar angles (0, <f>) of a real quasi-parallel interplanetary shock. The location 
of the solution of magnetic coplanarity (MC) velocity coplanarity (VC) ‘+\ 
Abraham-Shrauner (AS) ‘A’, Lepping-Argentiero (LA) *□’ and our solution (VS) 
are indicated. 


50 



0 


180 


60 120 
^ (deg) from X 

Figure 5c. Uniqueness’ contour plots of the log x 2 (0, 0) function versus the shock normal 
polar angles (0, 0) of a real planetary bow shock. The location of the solu- 
tion of magnetic coplanarity (MC) **’, velocity coplanarity (VC) V, Abraham- 
Shrauner (AS) ‘A’, Lepping-Argentiero (LA) *□’ and our solution (VS) *•’ are 
indicated. 
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Figure 6. ‘Uniqueness’ plots of the log x 2 (p) function versus the normalized density p = 
p/p* of a real a) quasi-perpendicular and b) quasi-parallel interplanetary shocks 
and c) planetary bow shock. The normalization constant p* is the value of the 
asymptotic density obtained by the iterative gradient scheme at the minimum of 
of the x function for each side of the shock. Furthermore, the plots on panel 
b show peaks located at p = 3.9 and 2.3 for the upstream and downstream 
sides of the shock corresponding to the singularity p = 47 tG 2 /B 2 associated 
with tangential, contact or rotational discontinuities. 
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Table 2a. Results of the analysis of the November 27, 1977 
interplanetary shock.- 
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Table 2b. Results of the analysis of the January 29, 1978 

interplanetary reverse shock. 
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Table 2c. Results of the analysis of the planetary bowshock. 
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Table 2c) gives 5679.8 ev/cc which has about 15? deviation from the observed 
value. This discrepancy can be explained by considering the errors incurred 
in the evaluation of the predicted pressure, since its calculation depends 
mostly in the poorly determined asymptotic magnetofluid variables. A crude 
estimate of the error bounds in the pressure jump due to uncertainties in 
the asymptotic magnetofluid variables yield ±970 ev/cc. It is clear then, 
that the predicted pressure jimp encompasses within this uncertainty the 
observed pressure jimp across the shock. Scudder et al . [1985] have also 
reported the pressure jump using a somewhat larger data interval. Their 
results are consistent within the error bounds to the values reported in 
this paper. 


5. Summary and Conclusions 

We have presented and demonstrated the utility of a new, fast, iterative 
method to determine the geometrical characteristics of a shock using the 
plasma and magnetic field observations together with a subset of a 
Rankine-Hugoniot model equations. The method exploited a new vector space 
that is separable, and unlike other methods contains a smaller number of 
non-linear unknown variables. An important aspect of the procedure is that 
'uniqueness' (or lack thereof) of the solutions can be demonstrated by 
either analytical or by graphical methods. To the best of our knowledge, 
this is the first time that 'uniqueness' has been demonstrated for the shock 
geometry solution. In so doing we also have illustrated the possible ways 
in which higher order non-linear techniques can obtain a misleading 
solution . 

The analysis we have presented indicates that, unlike extant methods, 
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this new iterative scheme is reliable at all 6g n ~angles regardless of the 

shock strength, geometry and direction of propagation relative to the 

ambient flow. The results in Tables la-c and 2a~c for synthetic and real 

shocks respectively, demonstrate the reliability and accuracy of the method 

in comparison to other procedures. A virtue of this method which indicates 

the well-conditioning of the approach is the lack of singular behavior for 

the extreme situations such as the purely perpendicular (0 O = 90°) and 

m 

parallel (0g = 0°) shock. Our analysis also indicates that the 

uncertainties in each set of parameters in the least squares sequence is 
smaller for the shock normal polar angles (i. e , the shock normal n) and 
increases for the specification of the asymptotic magnetofluid variables. 
This implies that the determination of the asymptotic states is more 
sensitive to errors in the observations. On the other hand, techniques such 
as magnetic coplanarity, velocity coplanarity and the Abraham-Shrauner mixed 
data pre-averaged methods select a priori these states to determine the 
shock normal and in doing so their shock normal calculation will be equally 
affected by these uncertainties. 

The comparison of shock parameters as obtained by different techniques 
indicates that seme of the other methods are reliable for particular shock 
geometries. In the case of perpendicular shocks, Abraham-Shrauner (AS) and 
velocity coplanarity (VC) methods gives good results for the shock geometry. 
On the other hand, magnetic coplanarity (MC) cannot describe the shock 
geometry of perpendicular shocks since its expression is singular as 0g n 
approaches 90°. Similarly the Lepping and Argentiero method cannot 
reasonably converge for even quasi-perpendicular shocks because its solution 
depends on the nearly singular expression of magnetic coplanarity. For 
parallel shocks we find that neither the MC, the LA nor the AS methods can 
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determine an accurate shock geometry. Again, this is because these methods 
are singular as 6g n approaches 0°. Generally all the techniques give 
reasonably good results for oblique shocks except for the approximate VC 
method which was demonstrated to fail in this geometry when the flow 
velocity was not aligned with the shock normal vector. 

There still remains various aspects on the determination of the shock 
geometry which deserve some consideration, however they can be difficult to 
implement. From the point of view of non-linear optimization, it is 
possible to incorporate the expression of the scalar pressure jump 
condition, even in the absence of temperature measurements, into the least 
squares normal equation for the shock normal polar angles determination. 
This condition will act as a constraint or penalty function and its effect 
will be to eliminate some of the unphysical solutions of the problem. 
Unfortunately, the analytical representation of this penalty function is not 
clear . 

An important application that resulted from our solution is the 
determination of various frames of references, such as the deHofftaan-Teller 
frame (HTF) [deHoffman and Teller, 1950] and the normal incidence frame 
(NIF) since their calculation depends on the shock normal, speed, 
conservation constants and the asymptotic magnetofluid states [Scudder et 
al., 1985]. With the availability of a technique that determines the 
optimal conserved fluxes at the shock, there is now a viable way to estimate 
these quantities which heretofore were expressed as functionals of the 
poorly determined 0g n values. For example, the deHoffman-Teller 
transformation velocity can be written either as 

V = 4 *'" tane Bn1 2 
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or as 



in terms of the conserved quantities of higher quality than the state 
variables . 
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Appendix A 


The analysis of the non-linear system of equations, such as for instance 
the Rankine-Hugoniot conservation equations (8) -- (11), is conveniently 

accomplished by means of the generalized inverse method. The application of 
this method to non-linear systems has been previously discussed, e.g. 
Jackson [1972], Bard [1974] and Lanczos [1961]. The generalized inverse 
method is a matrix formulation of the least squares problem where the 
fundamental equation to be solved is represented as 


A Ap = AY ( A 1 ) 

where AY :: Y - f(x^; p^° ) is a vector of length N* (i. e. i= 1, N') 

representing the difference between the observations Y and the model 

prediction and A is a matrix N' x M formed by the partial derivatives (i. e. 

the Jacobian) of the model equations with respect to the model unknown 

parameters (i. e. j = 1, M) evaluated at the initial guess. 

The solution of the normal equation (A1) is equivalent to the least 

2 ■> 

squares minimization method of the objective x (p) (i.e. the chi-square) 

function. This function is generally defined as 


9 N' 1 

X 2 (p) = 2 — 


2 (? i " VV p j )} 


(A2) 


i = 1 a 

where a represents the standard deviation of the observations. Equation 

(A2) gives a measure of how well the model equations represented by ?(x^; 

p.) fits the observations indicated by the vector Y. . In the matrix 
J 

2 

formalism, the minimization of the x function is analogous to the 
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determination of the optimun parameters that minimize the function 

X 2 (p) = r T r (A3) 

where r is the residual vector given by r = (Mp - A^j and r is the 

2 

transpose vector. Generally, the objective function x is normalized by the 
number of degrees of freedom v of the system. The number of degrees of 
freedom is defined as the difference between the total number of data points 
N' and the number of unknown parameters per model equation (M/L) (i. e. v = 
N' - M/L). Since the minimization of equation (A3) and the generalized 
method solution of equation (A1) have been shown to be mathematically 
equivalent [Lanczos, 1961; Jackson, 1972; Bard, 1974] we shall instead 
proceed with the application of the later method to the linearized matrix 
equation (A1). The reader is refered to the mentioned papers (and 
references therein) for the theoretical aspects of these methods. 

The matrix formulation of the generalized inverse method utilizes the 
singular value decomposition of Lanczos [Lanczos, 1961; Jackson, 1972]. 
This approach requires the estimation of the eigenvalues and eigenvectors 
associated with the matrix A in (A1). This approach is convenient when the 
matrix A is well conditioned in the sense that its eigenvalues are large and 
the iteration scheme will require short steps in the parameter space, 
keeping the linearization well inside its region of validity. However, if 
the matrix A is close to being ill-conditioned, which implies that some of 
its eigenvalues are zero or numerically very small, the solution vector Ap 
will take large steps in the parameter space that may well be outside of the 
region where the linearization is appropriate. This iterative process may 
then diverge unless some method of limiting the iterative step size is 
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employed. Two generally recognized options are used in this case. One 

option requires constructing a solution from the contribution of only the 

larger eigenvalues as suggested by Lanczos [1961] and Jackson [1972]. 

Although this procedure is reasonably appropriate, it requires the 

monitoring of the eigenvalues at each step in the iteration process making 

it slow. A second option, that we consider more practical and that can be 

easily implemented is to follow the technique known as the 

Marquardt-Levenberg* s algorithm [Levenberg, 1944; Marquardt, 1963; Bard, 

1974; Lawson and Hanson, 1974]. With this method the stability of the 

iterative procedure is improved by limiting the step size (more sensitive in 

the direction corresponding to the small eigenvalues) by introducing what is 

2 

known as a "cut-off" eigenvalue or Marquardt parameter a . Furthermore, with 
this "cut-off" eigenvalue, fast and accurate convergence is invoked and the 
need to monitor the small eigenvalues at each step of the iteration is 
avoided . 

The solution, then, of equation (A1) is now given by 
+ <*.1 -> 

Ap = A AY (A4) 

*” s 

where A is the generalized inverse defined by 

“*o 

Ag" 1 = (H 4- a 2 B)” 1 A T (A5) 

T T 

where A is the transpose matrix, H = A A is the approximate Hessian matrix 

which is positive definite and of size M x M, B is a diagonal matrix whose 

elements coincide with the diagonal elements of H if I~L .. £ 0 and with the 

2 

unit matrix I if = 0. The parameter « is the Marquardt parameter and 
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its size controls not only the step size but also the contribution of the 

small eigenvalues to the solution at each iteration step. 

■+ *► . 

In general the quantities in the vector Fyx^ py represent entities 
having different physical dimensions. For example, in the shock normal case 
in section 3a the function py is a vector of seven components 

representing the normal component of the magnetic field, the components of 
the tangential momentun flux and electric field in an arbitrary coordinate 
system. Since these quantities are constructed from the magnetic field and 
plasma observations, it is clear that some of the observations may be known 
to be less reliable than others, and we want to be certain that our 
parameter estimates will be less influenced by those than by the more 
accurate ones. For this reason it is convenient to weight equation (A1) 
before the parameters are estimated. After all, we cannot escape from the 
statistical nature of the observed data. One way of weighting the system of 
equations (A1) is by constructing the standard deviations associated with 
the physical variables of the Rankine-Hugoniot system. If the observations 
are statistically independent we define a diagonal matrix W = ( 1 /a^ ) of 
size N' x N' from the standard deviations. Operating on the normal equation 
( A 1 ) we have the solution 


AP* = (A 1 W T W A + a 2 B) _1 A T W T W A? 


(A6) 


2 2 -> y “p 

At this point the x function can be generalized to be x = r W W r, 


Let us now address the problem of the reliability and precision of the 


model parameters. It is not enough to compute a vector solution p without 
a simultaneous estimate of the error bounds in the parameters determined. 
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One way of expressing the reliability of the solution is by constructing 
what is called the resolution matrix [Lanczos, 1961; Jackson, 1972] given by 



The degree to which the R matrix approximates the identity matrix is a 
measure of the resolution obtainable from the data for each parameter. If 
the matrix R is nearly diagonal, then each parameter is a weighted sum of 
the others . 

To estimate the error bounds on the obtained parameters p , it is 

necessary to assume a statistical uncertainty distribution for them. This 

kind of test are exact only if the measurement errors do indeed follow such 

a distribution. Since in general such a distribution is unknown, a more 

practical way of obtaining the error bounds in the parameter space is to 

2 

consider the departure of the objective (risk) function x (p) from the 

2 -> 

obtained optimal value x ^P ) [Bard, 1974] as follows 


|x 2 (p) - X 2 (P*)| < e (A7) 

where e is the largest difference that one is willing to consider 
insignificant (i. e. the indifference region). Therefore we have no reason 
to prefer p over any other value of p for which (A7) is satisfied. The 
region enclosed in (A7) is named the indifference region. In a small 
neighborhood of p we can now approximate x (p) by a Taylor series expansion 


x 2 (p) 


X (P ) 


where 6p = p - p 


+#t ->■ i +t ® + 

+ q <$P + -y «P H 8p (A8) 


and q and H are the gradient vector and the Hessian 
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matrix of the x function respectively, evaluated at p = p . If p is an 

2 -*• -*.# 

optimal extremum of x (p). then q must vanish. We can now answer the 
question of the error bounds in the parameter p because equation (A8) 
properly written represents an M-dimensional ellipsoid whose principal axes 
(or eigenvalues of H) are a measure of these errors. Note that equation 
(A8) can now be written as 

<5p T H 6 P < 2e (A9) 

This is easily seen by noticing that equation (A9) can be formulated as an 
eigenvalue problem of the form 

H 6p = A 6p 


* 

where A ia a diagonal matrix of the eigenvalues of H . Thus operating by 
6p at both sides of this equation and using (A9) we get 

"T 

6P A Sp = 2e (A10) 

Equation (A10) states that the length of each vector component of 6p is 
proportional to / (2e/X) where X is its corresponding eigenvalue and the 
eigenvector represents the principal axis of the multidimensional ellipsoid. 
The largest axis (smallest eigenvalue) defines the worst-determined 
direction in p space and the shortest axis (largest eigenvalue) defines the 
best determined direction. Thus the solution of (A10) gives a reliable 
measure of the errors in the parameter determined. 

In a different way, if one does not have a good measure of the 
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indifference region e, it is possible to adopt an * ad hoc* error 
distribution, such as for example, the normal distribution and determine the 
confidence region e for p provided that the covariance matrix C of the 
errors of the observations is known [Bard, 1974; Jackson, 1972; Scheffe, 
1959!S. 

Acknowledments . We thank H. S. Bridge, N. F, Ness, K. W. Ogilvie and C 
Russell, principal investigators of the Voyager and ISEE--1 plasma and 
magnetic field experiments, respectively, for permission to use the plasma 
and magnetic field data. We would especially like to express our 

appreciation to R. Lepping and M. Acuna for the many stimulating 
discussions, their support and for allowing us to make use of their 
algorithm. We are greatful to E. Sittler for supplying the Voyager electron 
information and to Y. C. Whang for allowing us to use his shock algorithm. 
We also thank L. Burlaga, A. Lazarus and M. Goldstein for their continuous 
support . 


66 



References 


Abraham-Shrauner , B., Determination of magnetohydrodynamic shock normals, J » 
Geophys . Res ., 77, 736, 1972. 

Abraham-Shrauner , B. and S, H. Yun , Interplanetary shocks seen by Ames 
plasma probe on Pioneer 6 and 7. J . Geophys. Res . , 8 1 , 2097, 1976. 

Acuna, M. H. and R. P. Lepping, Modification to shock fitting program, J, 
Geophys. Res. , 89, 11004, 1984. 

Akhiezer, A. I., I. A. Akhiezer, R. V. Polovin, A. G. Sitenko and K. N. 
Stepanov, Plasma Electrodynamics; Linear Theory , Vol. I, Pergamon Press, 
New York, 1975. 

Armstrong, T. P. , M. E. Pesses and R. B. Decker, Shock drift acceleration, 
in Collisionless Shock Waves in the He liosphere , edited by E. C. Stone 
and B. Tsurutani, AGU Geomonograph, in press, 1985. 

Bard, Y. , Nonlin ear Parameter Estimation , Academic Press, New York, 1974. 

Boyd, T. J. and J. J. Sanderson, Plasma Dynamics , Barnes and Noble Inc., New 
York, 1969. 

Burlaga, L. F, , Hydromagnetic waves and discontinuities in the solar wind, 
Spac. Sci. Rev . , J_2, 600, 1971. 

Chao, J. K. , Interplanetary Collisionless Shock Waves, Rep. CSR TR-70-3, 
Massachusetts Institute of Technology, Cambridge, Mass., 1970. 

Chao, J. K. and K. C. Hsieh, On determining magnetohydrodynamic shock 
parameters © Bn and M A , Planet, Space Sci. , 32 , 641, 1984. 

Colburn, D. C. and C. P. Sonett, Discontinuities in the solar wind, Space 
Sci . Rev , , 5 , 439, 1966. 

deHoffman, F. and E. Teller, Magnetohydrodynamic shocks, Phys, Rev. , 80, 
692, 1950. 


67 



Edminston, J. P. and C. F. Kennel, A parametric survey of the first critical 
Mach number for fast MHD shock, J. Plasma Phys. , 32 , 429, 1985. 

Forman, M. A. and G. M. Webb, Acceleration of energetic particles, in 
Collisionless Shock Waves in the Heliosphere , edited by E. C. Stone and 
B. Tsurutani, AGU Geomonograph, in press, 1985. 

Goodrich, C. C. and J. D. Scudder, The adiabatic energy change of plasma 
electrons and the frame dependence of the cross-shock potential at 
collisionless magnetosonic shock waves, J. Geophys. Res ., 89 , 6654, 

1984. 

Gosling, J. T. , M. F. Thomsen, S. J. Bame, W. C. Feldman, G. Paschmann and 
N. Sckopke, Evidence for specularly reflected ions upstream from 
quasi-parallel bow shock, Geophys. Res. Lett ., 9 _, 1333, 1982. 

Greenstadt, E., V. Formisano, C. Goodrich, J. T. Gosling, M. Lee, M. Leroy, 
M. Mellot, K. Quest, A. E. Robson, P. Rodriguez, J. Scudder, J. Slavin, 
M. Thomsen, D. Winske and C. S. Wu , Collisionless Shock Waves in the 
Solar Terrestrial Enviroment, Solar Terrestrial Physics; Present and 
Future , edited by D. M. Butler and K. Papadopoulos , NASA publication 
1120, 1984. 

Jackson, D. D. , Interpretation of Innacurate, Insufficient and Inconsistent 


data, Geophys. J. R. 

Astr . Soc . , 

28, 

97, 

1972. 

Kennel , C. F. , J. P. 

Edmin ston 

and 

T. 

Hada, A quarter century of 


collisionless shock research, in Collisionless Shock Waves in the 
Heliosphere , edited by E. C. Stone and B. Tsurutani, AGU Geomonograph, 
in press, 1985. 

Lanczos, C. , Linear Differential Operators , D. Van Nostrand, Princeton, N. 
J., 1961. 

Landau, L. D. and E. M. Lifshitz, Electrodynamic of Continuous Media, 


68 



Pergaraon Press, pp 224-233, New York, I960. 


Lawson, C. L. and R. J. Hanson, Solving Least Squares Problems , Prentice 
Hall Inc, New Jersey, 197*1. 

Lepping, R. P. and P. D. Argentiero, Single spacecraft method of estimating 
shock normals, J. Geophys. Res ., 76 , 4349, 1971* 

Lepping, R. P. , Influence of thermal anisotropy on best-fit estimates of 
shock normals, J, Geophys. Res ., 77 , 2957. 1972. 

Levenberg, K. , A method for the solution of certain nonlinear problems in 
least squares, Qua rt. Appl. Math ., 2 , 164, 1944. 

Marquardt, D. W. , An algorithm for least squares estimation of nonlinear 
parameters, SIAM J . , 1J_, 431, 1963. 

Ogilvie, K. W. and L. F. Burlaga, Hydromagnetic shocks in the solar wind. 
Solar Phys. , 8, 422, 1969. 

Russell, C. T. , M. M. Mellott, E. J. Smith and J. H. King, Multiple 

spacecraft observations of interplanetary shocks: Four spacecraft 
determination of shock normals, J. Geophys. Res . , 88, 4739,1983a. 

Russell, C. T., J. T. Gosling, R. D. Zwickl and E, J. Smith, Multiple 

spacecraft observations of interplanetary shocks: ISEE three dimensional 
plasma measurements , J. Geophys. Res. , 88 , 9941, 1983b. 

Scheffe, H. , The Analysis of Variance , John Wiley & Sons Inc., New York, 
1959 . 

Scudder , J. D. , L. F, Burlaga and E. W. Greenstadt, Scale lengths in 

quasi-parallel shocks, J. Geophys. Res ., 89, 7545, 1984. 

Scudder, J. D. , C. Lacombe, A. Mageney, C. C. Harvey, J. T. Gosling, G. 
Paschmann, C. T. Russell, T. L. Aggson and R. Anderson, The resolved 

layer of a collisionless , high g, supercritical, quasiperpendicular 

shock wave." geometry, current system and scales, to be submitted to Jj, 


69 



Geophys. Res., 1985. 


Sonnerup, B. U. 0., Acceleration of particles reflected at a shock front, J. 
Geophys. Res. , 74 , 1301, 1969. 

Tidman, D. A. and N. A. Krall, Shock Waves in Collisionless Plasmas , John 
Wiley and Sons Inc., New York, 1971. 

Vinas, A. F. , M. L. Goldstein and M. H. Acuna, Spectral amalysis of 
magnetohydrodynamic fluctuations near interplanetary shocks, J. Geophys. 
Res. , 89, 3762, 1984. 

Whang, Y. C., F. Fei and H. Du, Transport of interplanetary fluctuations to 
the magnetopause, submitted to J. Geophys. Res ., 1985. 

Wu , C. S. , A fast Fermi process: Energetic electrons accelerated by a nearly 
perpendicular bow shock, J. Geophys. Res., 89 . 8857, 1984. 


A. F. Vinas and J. D. Scudder, Code 692, Laboratory for Extraterrestrial 
Physics, NASA/Goddard Space Flight Center, Greenbelt, MD, 20771* 


70 



1. Report No. 2. Government Accession No. 3. Recipient's Catalog No. 

TM-86214 


4. Title and Subtitle 5. Report Date 

FAST AND OPTIMAL SOLUTION TO THE ‘RANKINE- May 1985 

HUGONIOT PROBLEM’ 6. Performing Organization Code 


/. Author(s) 

Adolfo F. Vinas, Jack D. Scudder 


9. Performing Organization Name and Address 
NASA/Goddard Space Flight Center 
Greenbelt, Maryland 20771 


8. Performing Organization Report No. 

85B0392 


10. Work Unit No. 


1 1 . Contract or Grant No. 


12. Sponsoring Agency Name and Address 



tary Notes 


13. Type of Report and Period Covered 


Technical Memorandum 


14. Sponsoring Agency Code 


16. Abstract 

A new, definitive, reliable and fast iterative method is described for determining the geo- 
metrical properties of a shock (i.e. 0 Bn , ~n, V g and M A ), the conservation constants and the 
self-consistent asymptotic magnetofluid variables, that uses the three dimensional magnetic 
field and plasma observations. The method is well conditioned and reliable at all 0 Bn angle; 
regardless of the shock strength or geometry. Explicit proof of ‘uniqueness’ of the shock 
geometry solution by either analytical or graphical methods is given. The method is applied 
to synthetic and real shocks, including a bow shock event and the results are then compare! 
with those determined by preaveraging methods and other iterative schemes. A complete 
analysis of the confidence region and error bounds of the solution is also presented. 


17. Key Words (Selected by Author(s)) 
Rankine-Hugoniot , interplanetary 
shocks, tangential, rotational, 
contact discontinuity 


18. Distribution Statement 
Unclassified-Unlimited 
STAR: #75 and #88 


Unclassified 


"For sale by the National Technical Information Service, Springfield, Virginia 


20. Security Classif. (of this page) 

21 . No. of Pages 

Unclassified 

73 



GSFC 25-44 (10/7/) 






















